Análisis de Datos I
Tarea 5
Ejercicio 1.
Para la tabla de datos eurodist que viene con el
paquete datasets el cual contiene las distancias entre algunas de las
ciudades más importantes de Europa ejecute un Escalamiento
Multidimensional y luego compare el resultado con el mapa de
Europa.
Primero se obtienen los datos de las distancias apartir del paquete
datasets.
Una vez que se cuenta con las distancias, se ejecuta el ACP o Multidimensional Scalling, a partir de la matriz de distancias.
## $points
## [,1] [,2]
## Athens 2290.274680 1798.80293
## Barcelona -825.382790 546.81148
## Brussels 59.183341 -367.08135
## Calais -82.845973 -429.91466
## Cherbourg -352.499435 -290.90843
## Cologne 293.689633 -405.31194
## Copenhagen 681.931545 -1108.64478
## Geneva -9.423364 240.40600
## Gibraltar -2048.449113 642.45854
## Hamburg 561.108970 -773.36929
## Hook of Holland 164.921799 -549.36704
## Lisbon -1935.040811 49.12514
## Lyons -226.423236 187.08779
## Madrid -1423.353697 305.87513
## Marseilles -299.498710 388.80726
## Milan 260.878046 416.67381
## Munich 587.675679 81.18224
## Paris -156.836257 -211.13911
## Rome 709.413282 1109.36665
## Stockholm 839.445911 -1836.79055
## Vienna 911.230500 205.93020
##
## $eig
## [1] 1.953838e+07 1.185656e+07 1.528844e+06 1.118742e+06 7.893472e+05
## [6] 5.816552e+05 2.623192e+05 1.925976e+05 1.450845e+05 1.079673e+05
## [11] 5.139484e+04 -3.259629e-09 -9.496124e+03 -5.305820e+04 -1.322166e+05
## [16] -2.573360e+05 -3.326719e+05 -5.162523e+05 -9.191491e+05 -1.006504e+06
## [21] -2.251844e+06
##
## $x
## NULL
##
## $ac
## [1] 0
##
## $GOF
## [1] 0.7537543 0.8679134
Finalmente se gráfica a los individuos y se muestra una comparación con el mapa verdadero.
Ubicaciones apartir del MDS:
x <- MDS$points[,1]
y <- -1*MDS$points[,2]
plot(x, y, xlab="Componente 1", ylab="Componente 2",main="MDS",pch = 19)
text(x, y, labels = labels(matriz_distancias), cex=0.85 , pos = 1) Mapa real:
Como es posible observar, los resultados del MDS son muy similares al mapa real, con tan solo un para de diferencias que son atribuibles a no tener las distancias exactas o por lo contrario, al mapa con el que se cuenta.
Ejercicio 2.
Programe el algoritmo para Escalamiento Multidimensional
(Mutidimensional Scaling) visto en clase y con la tabla de estudiantes
compare con cmdscale() de R y veri que los resultados con
un ejemplo.
A continuación, se presenta una función que realiza el MDS paso a paso.
MDS_algoritmo <- function(X, distancias) {
# Paso 1: Con las distancias, se obtiene la matriz B
n <- nrow(X) # cantidad de individuos
B_matriz <- matrix(NA, n, n)
suma3 <- sum(distancias^2)
for (r in 1: n) {
suma1 <- 0
suma2 <- 0
for (s in 1: n ) {
suma1 <- sum((distancias[,s])^2)
suma2 <- sum((distancias[r,])^2)
B_matriz[r,s] <- -(1/2)*(distancias[r,s]^2 - (1/n)*(suma1 + suma2)
+ (1/n^2)*suma3)
}
}
# Paso 2 : Se calculan los valores y vectores propios de B.
B_e <- eigen(B_matriz)
valores_propios <- B_e$values
# Los valores propios muy pequeños cercanos a cero, se ignoran para el cálculo
# de las coordenadas. Se define un umbral para los valores propios
umbral <- 1e-10
# Se filtran los valores propios mayores que el umbral
valores_propios_filtrados <- valores_propios[valores_propios > umbral]
# Se obtienen los vectores propios asociados a los valores propios filtrados
m <- length(valores_propios_filtrados)
vectores_propios <- B_e$vectors[, 1:m]
# Paso 3 : Se calculan las coordenadas
coordenadas <- matrix(NA, n, m)
rownames(coordenadas) <- rownames(X)
for(i in 1: n) {
for(j in 1: m){
coordenadas[i,j] <- (sqrt(valores_propios_filtrados[j])*
vectores_propios[i,j])
}
}
# Paso 4: Graficar
x <- coordenadas[,1]
y <- coordenadas[,2]
grafico <- plot(x, y, xlab="Componente 1", ylab="Componente 2",
main="MDS algoritmo",pch = 19)
text(x, y, labels = row.names(X), cex = 0.85 , pos = 1)
resultado <- list("points" = coordenadas, "eig" = valores_propios,
"plot"= grafico)
return(resultado)
}Comparación con cmdscales()
Primeramente, se carga la base de datos de estudiantes.
Se obtiene las distancias de cada par de individuos.
Se verifican los resultados obtenidos por ambas funciones
## $points
## [,1] [,2] [,3] [,4] [,5]
## Lucia 0.76471745 1.5817637 1.11186219 0.039908252 -0.0007926494
## Pedro -1.66887794 -1.3919656 0.09067929 -0.053941739 0.1128975404
## Ines -1.57822841 -0.2994960 0.48752985 0.552652786 -0.1377925976
## Luis 2.60701317 -1.3202040 -0.46230941 0.784004734 0.0524431575
## Andres 1.43877557 1.3356687 -0.67985389 -0.189277341 -0.1117444886
## Ana -2.34790534 -0.3880845 -0.12895699 0.007628838 -0.0253542859
## Carlos 0.89372557 1.5189012 -0.38893244 -0.124247926 -0.0093618642
## Jose -2.64984571 -0.4254636 -0.46447580 -0.316066094 -0.0162142654
## Sonia 2.62959083 -2.1833951 0.40705140 -0.658738298 -0.0317927290
## Maria -0.08896518 1.5722752 0.02740580 -0.041923214 0.1677121823
##
## $eig
## [1] 3.498309e+01 1.793415e+01 2.708155e+00 1.511504e+00 7.710194e-02
## [6] 3.775779e-15 2.673082e-15 2.327544e-15 -9.116253e-16 -3.093847e-15
##
## $x
## NULL
##
## $ac
## [1] 0
##
## $GOF
## [1] 1 1
x <- MDS_estudiantes$points[,1]
y <- MDS_estudiantes$points[,2]
plot(x, y, xlab="Componente 1", ylab="Componente 2",
main="MDS",pch = 19)
text(x, y, labels = row.names(estudiantes_datos), cex = 0.85 , pos = 1)Las coordenadas de los individuos son iguales. Por otro lado, los primeros 5 valores propios dados por ambos son los mismos. En cuanto a los restantes 4, estos difieren,pues, al ser números cercanos a cero por cuestiones de cálculo de la función cmdscale la precisión númerica cambia en comparación con eigen.Por último, el plano principal es el mismo.
Ejercicio 3.
Pruebe la ecuación (1.3) de la filminas vistas en clase, la
cual establece una fórmula para el cálculo de brs.
Ejercicio 4.
Pruebe el Teorema que establece la relación entre los valores
y vectores propios del algoritmo anterior en relación a los valores y
vectores propios del ACP.
Ejercicio 5.
Del documento adjunto MDS no métrico.pdf leer la sección MDS no métrico y replicar en R los dos ejemplos que ah se presentan.
Ejemplo 1.
nombres_fila_columna1 <- c("M", "B", "V", "S", "SS", "LC")
m <- matrix(c(
0, 627, 351, 550, 488, 603,
627, 0, 361, 1043,565, 1113,
351, 361, 0, 567, 564, 954,
550, 1043,567, 0, 971, 950,
488, 565, 564, 971, 0, 713,
603, 1113,954, 950, 713, 0),
nrow = 6, byrow = TRUE)
m <- as.data.frame(m)
rownames(m) <- nombres_fila_columna1
colnames(m) <- nombres_fila_columna1
m <- as.matrix(m)Calculamos Q
Q <- matrix(0, nrow = 6, ncol = 6)
sum1 <- c()
i <- 0
for (s in 1:6) {
for (r in 1:6) {
i <- i + (m[r,s])^2
sum1[s] <- i
}
i <- 0
}
sum2 <- c()
j <- 0
for (r in 1:6) {
for (s in 1:6) {
j <- j + (m[r,s])^2
sum2[r] <- j
}
j <- 0
}
sum3 <- 0
for (r in 1:6) {
for (s in 1:6) {
sum3 <- sum3 + (m[r,s])^2
}
}
for (r in 1:6) {
for (s in 1:6) {
Q[r,s] <- -(1/2)*(((m[r,s])^2) - (1/6)*sum1[s] - (1/6)*sum2[r] +
(1/6^2)*sum3)
}
}
Q ## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 11759.44 -39079.22 -17954.39 38559.11 -31804.89 38519.94
## [2,] -39079.22 303211.11 124211.44 -208389.56 73380.44 -253334.22
## [3,] -17954.39 124211.44 75532.78 60951.28 -39894.22 -202846.89
## [4,] 38559.11 -208389.56 60951.28 367858.78 -206103.72 -52875.89
## [5,] -31804.89 73380.44 -39894.22 -206103.72 162774.78 41647.61
## [6,] 38519.94 -253334.22 -202846.89 -52875.89 41647.61 428889.44
Encontramos los valores y vectores propios de Q
propios <- eigen(Q)
valores_propios <- propios$values
valores_propios <- valores_propios[-c(5,6)]
vectores_propios <- propios$vectors
vectores_propios <- vectores_propios[,-c(5,6)]
valores_propios## [1] 737916.932 591058.808 59468.730 1035.536
## [,1] [,2] [,3] [,4]
## [1,] 0.0959729 -0.04429354 0.256942195 0.85659994
## [2,] -0.6270080 0.14004992 0.415521275 -0.15933263
## [3,] -0.2832133 -0.25835017 0.009409856 -0.31296904
## [4,] 0.2934011 -0.72163136 -0.220480281 -0.12846516
## [5,] -0.1240990 0.44165813 -0.781233037 0.08849908
## [6,] 0.6449462 0.44256702 0.319839993 -0.34433219
Calculamos las coordenadas
coordenadas <- matrix(0, nrow = 6, ncol = 4)
for (i in 1:6) {
for (j in 1:4) {
coordenadas[i,j] <- sqrt(valores_propios[j])*vectores_propios[i,j]
}
}
coordenadas[,1:2]## [,1] [,2]
## [1,] 82.44273 -34.05303
## [2,] -538.61294 107.67087
## [3,] -243.28619 -198.62051
## [4,] 252.03769 -554.79271
## [5,] -106.60360 339.54831
## [6,] 554.02232 340.24707
Graficamos
plot(coordenadas[,1], coordenadas[,2], col = "orange", pch = 16,
xlab = "Coordenada X", ylab = "Coordenada Y", main = "Plot de Coordenadas")
text(coordenadas[,1], coordenadas[,2],
labels = rownames(m), pos = 4, cex = 0.7, col = "black")Ejemplo 2.
nombres_fila_columna2 <- c("Atlanta", "Chicago", "Denver", "Houston",
"L. Angeles", "Miami", "N York",
"S Francisco", "Seattle", "Washington")
aire.dist <- matrix(c(0, 587, 1212, 701, 1936, 604, 748, 2139, 218, 543,
587, 0, 920, 940, 1745, 1188, 713, 1858, 1737, 597,
1212, 920, 0, 879, 831, 1726, 1631, 949, 1021, 1494,
701, 940, 879, 0, 1374, 968, 1420, 1645, 1891, 1220,
1936, 1745, 831, 1374, 0, 2339, 2451, 347, 959, 2300,
604, 1188, 1726, 968, 2339, 0, 1092, 2594, 2734, 923,
748, 713, 1631, 1420, 2451, 1092, 0, 2571, 2408, 205,
2139, 1858, 949, 1645, 347, 2594, 2571, 0, 678, 2442,
218, 1737, 1021, 1891, 959, 2734, 2408, 678, 0, 2329,
543, 597, 1494, 1220, 2300, 923, 205, 2442, 2329, 0),
nrow = 10, byrow = TRUE)
aire.dist <- as.data.frame(aire.dist)
rownames(aire.dist) <- nombres_fila_columna2
colnames(aire.dist) <- nombres_fila_columna2Se efectua un analisis clasico MDS metrico
## Warning in cmdscale(as.matrix(aire.dist), k = 9, eig = T): only 5 of the first
## 9 eigenvalues are > 0
## $points
## [,1] [,2] [,3] [,4] [,5]
## Atlanta -434.7588 724.22221 440.92520 0.18579147 -0.0125796337
## Chicago -412.6102 55.04016 -370.93031 4.39607623 12.6754984940
## Denver 468.1952 -180.65789 -213.57331 30.40856601 -9.5854646599
## Houston -175.5816 -515.22265 362.83981 9.48713270 -4.8603541375
## L. Angeles 1206.6772 -465.63705 56.52609 1.34143943 6.8086243199
## Miami -1161.6875 -477.98261 479.59934 -13.79783476 2.2781800452
## N York -1115.5609 199.79247 -429.66594 -29.39692855 -7.1366804644
## S Francisco 1422.6887 -308.65595 -205.51788 -26.06309834 -1.9830609621
## Seattle 1221.5351 887.20174 170.44892 -0.06998691 -0.0000894326
## Washington -1018.8972 81.89956 -290.65190 23.50884273 1.8159264311
##
## $eig
## [1] 9.213705e+06 2.199924e+06 1.082863e+06 3.322361e+03 3.858824e+02
## [6] -5.234966e-09 -9.323115e+01 -2.168535e+03 -9.090644e+03 -1.722963e+06
##
## $x
## NULL
##
## $ac
## [1] 0
##
## $GOF
## [1] 0.8781612 1.0000000
Se Calculan los autovalores
## [1] 9.213705e+06 2.199924e+06 1.082863e+06 3.322361e+03 3.858824e+02
## [6] -5.234966e-09 -9.323115e+01 -2.168535e+03 -9.090644e+03 -1.722963e+06
Se normalizan los dos primeros autovalores
## [1] 0.8018277
## [1] 0.9558842
Se muestran las coordenadas de las ciudades en las dos dimensiones
## [,1] [,2]
## Atlanta -434.7588 724.22221
## Chicago -412.6102 55.04016
## Denver 468.1952 -180.65789
## Houston -175.5816 -515.22265
## L. Angeles 1206.6772 -465.63705
## Miami -1161.6875 -477.98261
## N York -1115.5609 199.79247
## S Francisco 1422.6887 -308.65595
## Seattle 1221.5351 887.20174
## Washington -1018.8972 81.89956
Se dibujan las coordenadas de las ciudades en las dos dimensiones
par(pty="s")
plot(-aire.mds$points[,1],aire.mds$points[,2],type="n",xlab="Coordenada 1",
ylab="Coordenada 2",xlim = c(-2000,1500),ylim=c(-2000,1500))
text(-aire.mds$points[,1],aire.mds$points[,2],labels=row.names(aire.dist))otra manera de dibujar las coordenadas
plot(-aire.mds$points[,1], aire.mds$points[,2], col = "orange", pch = 16,
xlab = "Coordenada X", ylab = "Coordenada Y", main = "Plot de Coordenadas")
text(-aire.mds$points[,1], aire.mds$points[,2],
labels = rownames(aire.mds$points), pos = 3, cex = 0.6, col = "black")El fichero original de datos viene recogido ya en MASS:
## Fertility Agriculture Examination Education Catholic
## Courtelary 80.2 17.0 15 12 9.96
## Delemont 83.1 45.1 6 9 84.84
## Franches-Mnt 92.5 39.7 5 5 93.40
## Moutier 85.8 36.5 12 7 33.77
## Neuveville 76.9 43.5 17 15 5.16
## Porrentruy 76.1 35.3 9 7 90.57
## Broye 83.8 70.2 16 7 92.85
## Glane 92.4 67.8 14 8 97.16
## Gruyere 82.4 53.3 12 7 97.67
## Sarine 82.9 45.2 16 13 91.38
## Veveyse 87.1 64.5 14 6 98.61
## Aigle 64.1 62.0 21 12 8.52
## Aubonne 66.9 67.5 14 7 2.27
## Avenches 68.9 60.7 19 12 4.43
## Cossonay 61.7 69.3 22 5 2.82
## Echallens 68.3 72.6 18 2 24.20
## Grandson 71.7 34.0 17 8 3.30
## Lausanne 55.7 19.4 26 28 12.11
## La Vallee 54.3 15.2 31 20 2.15
## Lavaux 65.1 73.0 19 9 2.84
## Morges 65.5 59.8 22 10 5.23
## Moudon 65.0 55.1 14 3 4.52
## Nyone 56.6 50.9 22 12 15.14
## Orbe 57.4 54.1 20 6 4.20
## Oron 72.5 71.2 12 1 2.40
## Payerne 74.2 58.1 14 8 5.23
## Paysd'enhaut 72.0 63.5 6 3 2.56
## Rolle 60.5 60.8 16 10 7.72
## Vevey 58.3 26.8 25 19 18.46
## Yverdon 65.4 49.5 15 8 6.10
## Conthey 75.5 85.9 3 2 99.71
## Entremont 69.3 84.9 7 6 99.68
## Herens 77.3 89.7 5 2 100.00
## Martigwy 70.5 78.2 12 6 98.96
## Monthey 79.4 64.9 7 3 98.22
## St Maurice 65.0 75.9 9 9 99.06
## Sierre 92.2 84.6 3 3 99.46
## Sion 79.3 63.1 13 13 96.83
## Boudry 70.4 38.4 26 12 5.62
## La Chauxdfnd 65.7 7.7 29 11 13.79
## Le Locle 72.7 16.7 22 13 11.22
## Neuchatel 64.4 17.6 35 32 16.92
## Val de Ruz 77.6 37.6 15 7 4.97
## ValdeTravers 67.6 18.7 25 7 8.65
## V. De Geneve 35.0 1.2 37 53 42.34
## Rive Droite 44.7 46.6 16 29 50.43
## Rive Gauche 42.8 27.7 22 29 58.33
## Infant.Mortality
## Courtelary 22.2
## Delemont 22.2
## Franches-Mnt 20.2
## Moutier 20.3
## Neuveville 20.6
## Porrentruy 26.6
## Broye 23.6
## Glane 24.9
## Gruyere 21.0
## Sarine 24.4
## Veveyse 24.5
## Aigle 16.5
## Aubonne 19.1
## Avenches 22.7
## Cossonay 18.7
## Echallens 21.2
## Grandson 20.0
## Lausanne 20.2
## La Vallee 10.8
## Lavaux 20.0
## Morges 18.0
## Moudon 22.4
## Nyone 16.7
## Orbe 15.3
## Oron 21.0
## Payerne 23.8
## Paysd'enhaut 18.0
## Rolle 16.3
## Vevey 20.9
## Yverdon 22.5
## Conthey 15.1
## Entremont 19.8
## Herens 18.3
## Martigwy 19.4
## Monthey 20.2
## St Maurice 17.8
## Sierre 16.3
## Sion 18.1
## Boudry 20.3
## La Chauxdfnd 20.5
## Le Locle 18.9
## Neuchatel 23.0
## Val de Ruz 20.0
## ValdeTravers 19.5
## V. De Geneve 18.0
## Rive Droite 18.2
## Rive Gauche 19.3
## initial value 2.979731
## iter 5 value 2.431486
## iter 10 value 2.343353
## final value 2.338839
## converged
Las coordenadas son
## [,1] [,2]
## Courtelary 39.976648 -18.6652658
## Delemont -42.144880 -15.8354737
## Franches-Mnt -49.183472 -23.5017064
## Moutier 10.640047 -7.7948873
## Neuveville 35.713299 4.2241084
## Porrentruy -45.069549 -26.7640938
## Broye -55.273717 3.1907784
## Glane -58.641853 0.1769567
## Gruyere -55.274386 -11.8305230
## Sarine -46.800689 -18.2672079
## Veveyse -59.271433 -2.7666474
## Aigle 27.639864 18.9592510
## Aubonne 31.354629 26.3662771
## Avenches 32.353450 19.3520700
## Cossonay 31.342613 28.5036282
## Echallens 10.799996 26.1967959
## Grandson 40.595874 -1.7999750
## Lausanne 36.972991 -25.1351295
## La Vallee 50.122185 -23.5248585
## Lavaux 30.346198 30.8619325
## Morges 31.803525 18.6695352
## Moudon 33.621794 16.2266437
## Nyone 25.423297 7.4959696
## Orbe 34.176435 14.5941476
## Oron 30.076709 32.2309096
## Payerne 31.559866 17.6683208
## Paysd'enhaut 32.631374 27.0802839
## Rolle 28.769137 19.1245372
## Vevey 29.599603 -16.5009238
## Yverdon 33.286505 10.5935022
## Conthey -67.755640 18.1702425
## Entremont -66.402167 14.4997009
## Herens -68.766153 20.1878849
## Martigwy -63.370151 8.9517543
## Monthey -60.005854 -0.5353731
## St Maurice -62.872702 7.3903698
## Sierre -66.745414 16.5445662
## Sion -56.775840 -4.1215356
## Boudry 37.427364 -1.3796084
## La Chauxdfnd 41.684937 -30.7139972
## Le Locle 39.469769 -20.2862591
## Neuchatel 34.177795 -33.7884464
## Val de Ruz 37.403309 1.2494697
## ValdeTravers 42.204515 -17.0228055
## V. De Geneve 17.524061 -65.0238668
## Rive Droite -6.829233 -11.9126783
## Rive Gauche -7.514655 -31.3383738
Se dibujan los puntos
plot(swiss.mds$points, type="n", xlab="Coordenada 1", ylab="Coordenada 2")
text(swiss.mds$points, labels=as.character(1:nrow(swiss.x)))Otra manera de dibujar los puntos con los nombres de etiqueta
coordenadas <- swiss.mds$points
nombres <- rownames(swiss.mds$points)
plot(coordenadas[,1], coordenadas[,2], col = "orange", pch = 16,
xlab = "Coordenada X", ylab = "Coordenada Y", main = "Plot de Coordenadas")
text(coordenadas[,1], coordenadas[,2], labels = nombres, pos = 4, cex = 0.3,
col = "black")